# MODELO MULTIFASE – BUCARAMANGA
# (Clusterización + Localización + Evaluación de Cobertura)

cat("\014") # Borra la consola (en RStudio )
graphics.off() # Cierra todas las ventanas/dispositivos gráficos abiertos
# Antes de ejecutar el modelo, se eliminan todos los archivos
# generados en ejecuciones anteriores para evitar confusión
# entre resultados antiguos y nuevos.
OUT_DIR <- "outputs"

if (dir.exists(OUT_DIR)) {
  files_old <- list.files(OUT_DIR, full.names = TRUE)
  if (length(files_old) > 0) {
    file.remove(files_old)
  }
} else {
  dir.create(OUT_DIR)
}
# Verificación del entorno de ejecución:
#   Se asegura que el script se ejecute dentro del proyecto MODELO.Rproj,
#   de manera que las rutas relativas a las carpetas (DATA/, outputs/)
#   se resuelvan correctamente y se garantice la reproducibilidad del modelo.

if (!file.exists("MODELO.Rproj")) {
  stop("Este script debe ejecutarse abriendo 
       primero MODELO.Rproj (rutas relativas).")
}
# ------------------------------
# Paquetes y explicación de uso
# ------------------------------
# cluster     : daisy(Gower), pam, silhouette (clusterización sobre disimilitudes)
# factoextra  : fviz_dist, fviz_nbclust (codo/silhouette), fviz_dend (dendrograma)
# tidyr/dplyr : resúmenes, perfiles, transformaciones
# sf          : distancias métricas (cobertura), reproyección CRS
# ClustGeo    : soporte metodológico (integración distancia espacial + atributos)
# ade4        : is.euclid + cailliez (corrección a euclidianidad)
# fpc/scales/ggplot2/ape/units : utilidades y gráficos

ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
  invisible(sapply(pkg, require, character.only = TRUE))
}
packages <- c(
  "fpc","cluster","factoextra","tidyr","readxl","dplyr",
  "sf","ClustGeo","scales","ggplot2","ape","units","ade4","readr"
)
ipak(packages)
options(lifecycle_verbosity = "quiet") #suprimir mensajes informativos no críticos
# ------------------------------------------
# Lectura del ÚNICO input (CONSOLIDADO_BGA)
# ------------------------------------------
FILE_CONSOLIDADO <- file.path("DATA", "CONSOLIDADO_BGA.xlsx")
if (!file.exists(FILE_CONSOLIDADO)) {
  stop("No se encontró CONSOLIDADO_BGA.xlsx en la carpeta /DATA")
}

CONSOLIDADO_BGA <- readxl::read_excel(FILE_CONSOLIDADO)
message("[OK] CONSOLIDADO_BGA.xlsx cargado desde /DATA")
# --------------------
# Parámetros globales
# --------------------
# Parámetro de para el análisis de selección de kclústeres
K_MAX      <- 15 

# Parámetros para fusionar las matrices
W_ATR      <- 0.5 
W_ESP      <- 0.5 

# Sistemas de referencia usados para coordenadas y calcular distancias
CRS_GEO    <- 4326
CRS_PROJ   <- 3116

# Parámetros para estimar demanda potencial
WMIN <- 0.005  # 0.5%
WMAX <- 0.05   # 5%

# Parámetros para calcular el radio de cobertura
T_RTA_MIN <- 10
V_NOM_KMH <- 20

# Escenarios a revisar(selección de los kclústeres)
# Pre-estaablecido ya para los esquemas a revisar de acuerdo con el
# desarrollo del modelo con base en la sección 5.2.1 del libro: 
# Definición de los esquemas de localización
K_ELEGIDOS <- c(2,3,5)

# Directorio de salidas
OUT_DIR <- "outputs"
dir.create(OUT_DIR, showWarnings = FALSE)
# ---------------------------------------
# Definición de las funciones auxiliares
# ---------------------------------------
# Asigna etiquetas (nombre_barrios) a un objeto de tipo 'dist',
# preservando la coherencia entre la matriz y las observaciones
set_labels <- function(d, labels){
  stopifnot(inherits(d, "dist"), length(labels) == attr(d, "Size"))
  attr(d, "Labels") <- labels
  d
}
# Exporta un objeto de distancias ('dist') a un archivo CSV, convirtiéndolo
# previamente a matriz y redondeando los valores
write_dist_csv <- function(d, file, digits = 6){
  stopifnot(inherits(d, "dist"))
  M <- round(as.matrix(d), digits)
  utils::write.csv(M, file = file, row.names = TRUE)
}
# Calcula la moda de un vector numérico. En caso de existir múltiples modas,
# retorna todas concatenadas, separadas por '/'.
moda_numerica <- function(x){
  x <- x[!is.na(x)]
  if (!length(x)) return(NA_character_)
  tab   <- table(x)
  max_f <- max(tab)
  modos <- as.numeric(names(tab)[tab == max_f])
  paste(modos, collapse = "/")
}
# Calcula la moda del estrato socioeconómico de un conjunto de barrios.
# Se utiliza como valor representativo del clúster para el cálculo del radio de cobertura.
moda_estrato <- function(x){
  x <- x[!is.na(x)]
  if(length(x) == 0) return(NA_integer_)
  tb <- table(x)
  as.integer(names(tb)[which(tb == max(tb))][1])
}
# Realiza una normalización min–max sobre un vector numérico,
# retornando valores en el intervalo [0,1].
norm01 <- function(x){
  r <- range(x, na.rm = TRUE)
  if (r[1] == r[2]) return(rep(0, length(x)))
  (x - r[1]) / (r[2] - r[1])
}
# Normalización min–max usando límites externos (p.ej., rango global ciudad)
norm01_lim <- function(x, xmin, xmax){
  if (is.na(xmin) || is.na(xmax) || xmin == xmax) return(rep(0, length(x)))
  (x - xmin) / (xmax - xmin)
}
# Imprime objetos gráficos o resultados suprimiendo advertencias no críticas,
# con el fin de mantener la salida en consola limpia.
print_quiet <- function(p) suppressWarnings(print(p))
#------------------------------------------------------------
# Función auxiliar para la estimación de demanda potencial
  # Cálculo de demanda:
    #  - V_ESTRATO = (6 - estrato)/5
    #  - V_IPM = min–max(IPM)
    #  - V_DANO = min–max(%DAÑO)
    #  - V_BARRIO = promedio simple
    #  - Wi = WMIN + (WMAX-WMIN)*V_BARRIO
    #  - DEMANDA = (POB * Wi)-redondeado
#-------------------------------------------------------------  
calc_demanda <- function(df, col_pob, col_estrato, col_ipm, col_danio,
                         WMIN = 0.005, WMAX = 0.05){
  # Extracción y conversión de variables base
  P  <- as.numeric(df[[col_pob]])
  E  <- pmin(pmax(as.numeric(df[[col_estrato]]), 1), 6)
  IP <- as.numeric(df[[col_ipm]])
  D  <- as.numeric(df[[col_danio]])
  # Construcción de indicadores normalizados de vulnerabilidad
  V_ESTRATO <- (6 - E) / 5
  V_IPM     <- norm01(IP)
  V_DANO    <- norm01(D)
  # Índice de vulnerabilidad del barrio
  V_BARRIO  <- (V_ESTRATO + V_IPM + V_DANO) / 3
  V_BARRIO  <- pmin(pmax(V_BARRIO, 0), 1)
  # Factor de ponderación de la demanda y cálculo de demanda potencial
  Wi        <- WMIN + (WMAX - WMIN) * V_BARRIO
  DEMANDA   <- ceiling(P * Wi)
  # Retorno de resultados intermedios y demanda final
  data.frame(
    V_ESTRATO = V_ESTRATO,
    V_IPM     = V_IPM,
    V_DANO    = V_DANO,
    V_BARRIO  = V_BARRIO,
    Wi        = Wi,
    DEMANDA   = DEMANDA
  )
}
# --------------------------------------
# Entrada de datos y mapeo de columnas
# --------------------------------------
df0 <- as.data.frame(CONSOLIDADO_BGA)
df0$BARRIO <- as.character(df0$BARRIO)

# Asigna el nombre del barrio como identificador de fila del conjunto de datos
rownames(df0) <- as.character(df0$BARRIO) 
df0$BARRIO   <- NULL

# Correspondencia entre columnas del archivo de entrada y variables del modelo
map <- list(
  estrato   = "ESTRATO",
  ipm       = "IPM",
  poblacion = "POBLACIÓN",
  danio     = "% DAÑO",
  lon       = "X",
  lat       = "Y"
)

# Verifica que todas las columnas requeridas por el modelo estén presentes
# en el archivo de entrada antes de continuar con la ejecución
stopifnot(all(c(map$estrato, map$ipm, map$poblacion,
                map$danio, map$lon, map$lat) %in% names(df0)))

# -------------------------------
# Cálculo de demanda potencial 
# ------------------------------
dem <- calc_demanda(
  df         = df0,
  col_pob    = map$poblacion,
  col_estrato= map$estrato,
  col_ipm    = map$ipm,
  col_danio  = map$danio,
  WMIN = WMIN, WMAX = WMAX
)

df0 <- cbind(df0, dem)

barrios_dem <- data.frame(
  BARRIO  = rownames(df0),
  DEMANDA = as.numeric(df0$DEMANDA),
  X       = as.numeric(df0[[map$lon]]),
  Y       = as.numeric(df0[[map$lat]]),
  stringsAsFactors = FALSE
)
readr::write_csv(barrios_dem, file.path(OUT_DIR, "barrios_demanda_centroides_calculada.csv"))

message("[OK] DEMANDA calculada  y exportada: outputs/barrios_demanda_centroides_calculada.csv")
message(sprintf("[i] Parámetros DEMANDA: WMIN=%.3f%%, WMAX=%.3f%%", WMIN*100, WMAX*100))
# ------------------------------------------------
# Preparación y fusión de matrices de distancia
# ------------------------------------------------
#  - D0e: matriz de distancia multidemensional (Gower) corregida a euclidiana
#  - D1e: matriz de distancia espacial euclidiana (CRS métrico) 
#  - Dfused_e: matriz de distancia fusionada, utilizada como insumo del clustering PAM
preparar_y_fusionar <- function(df_base, w_espacial = W_ESP){
  
# (1) Preparación de variables para distancia socioeconómica
  # Conversión explícita de variables a formato numérico y
  # acotación del estrato a valores válidos [1,6].
  df <- df_base %>%
    dplyr::mutate(
      ipm_01        = as.numeric(.data[[map$ipm]]),
      danio_01      = as.numeric(.data[[map$danio]]),
      pob_raw       = as.numeric(.data[[map$poblacion]]),
      estrato_temp  = pmin(pmax(as.numeric(.data[[map$estrato]]), 1), 6)
    )
  # Definición del estrato como variable ordinal (requerido por Gower)
  df$estrato_ord <- ordered(df$estrato_temp, levels = 1:6)
  # Matriz de atributos utilizada para el cálculo de distancias socioeconómicas
  X <- df[, c("danio_01", "ipm_01", "pob_raw", "estrato_ord")]
  rownames(X) <- rownames(df_base)

# (2) D0
  # Se utiliza la métrica de Gower para manejar variables mixtas(numéricas y ordinales)
  D0  <- cluster::daisy(X, metric = "gower")
  D0d <- stats::as.dist(as.matrix(D0))
  D0d <- set_labels(D0d, rownames(X))
  
  # Verificación de euclidianidad y corrección mediante Cailliez si es necesario
  if (!suppressWarnings(ade4::is.euclid(D0d))) {
    D0e <- suppressWarnings(ade4::cailliez(D0d))
    c0  <- mean(as.numeric(D0e - D0d))
    message(sprintf("D0 (atributos) no euclidiana. Corrección Cailliez c* ≈ %.6f.", c0))
  } else {
    D0e <- D0d
    message("D0 (atributos) ya es euclidiana. No se aplica Cailliez.")
  }
  D0e <- set_labels(D0e, rownames(X))
  
# (3) D1
  # Extracción de coordenadas geográficas y transformación a
  # sistema de referencia proyectado (CRS 3116) para medir distancias métricas.
  df$`__lon` <- as.numeric(df_base[[map$lon]])
  df$`__lat` <- as.numeric(df_base[[map$lat]])
  
  pts      <- sf::st_as_sf(df, coords = c("__lon","__lat"),
                           crs = CRS_GEO, remove = FALSE)
  pts_proj <- sf::st_transform(pts, crs = CRS_PROJ)
  
  coords <- sf::st_coordinates(pts_proj)
  x <- coords[, "X"]
  y <- coords[, "Y"]
  
  # Cálculo explícito de distancias euclidianas(en kilómetros)
  DX <- abs(outer(x, x, "-"))
  DY <- abs(outer(y, y, "-"))
  M  <- sqrt(DX^2 + DY^2) / 1000
  diag(M) <- 0
  
  D1d <- stats::as.dist(M)
  D1d <- set_labels(D1d, rownames(X))
  message("Distancia espacial base calculada en CRS 3116 (euclidiana, km).")
  
  # Verificación de euclidianidad y corrección si aplica
  if (!suppressWarnings(ade4::is.euclid(D1d))) {
    D1e <- suppressWarnings(ade4::cailliez(D1d))
    c1  <- mean(as.numeric(D1e - D1d))
    message(sprintf("D1 (espacial) no euclidiana. Corrección Cailliez c* ≈ %.6f.", c1))
  } else {
    D1e <- D1d
    message("D1 (espacial) ya es euclidiana. No se aplica Cailliez.")
  }
  D1e <- set_labels(D1e, rownames(X))
  
# (4) Fusión de matrices de distancia
  # Combinación lineal controlada por el parámetro w_espacial.
  MF <- (1 - w_espacial) * as.matrix(D0e) +
    w_espacial      * as.matrix(D1e)
  diag(MF) <- 0
  
  Dfused <- stats::as.dist(MF)
  Dfused <- set_labels(Dfused, rownames(X))
  
  message(sprintf("Fusión de distancias: D = %.2f·D0e + %.2f·D1e",
                  (1 - w_espacial), w_espacial))
  
  # Verificación final de euclidianidad de la matriz fusionada
  if (!suppressWarnings(ade4::is.euclid(Dfused))) {
    Dfused_e <- suppressWarnings(ade4::cailliez(Dfused))
    cf <- mean(as.numeric(Dfused_e - Dfused))
    message(sprintf("Dfused no euclidiana. Corrección Cailliez final c* ≈ %.6f.", cf))
  } else {
    Dfused_e <- Dfused
    message("Dfused ya es euclidiana. No se aplica Cailliez.")
  }
  Dfused_e <- set_labels(Dfused_e, rownames(X))

# (5) Salida de la función
  # Retorna las matrices necesarias para las fases posteriores del modelo.
  list(
    Dfused_e = Dfused_e,  # matriz de distancia final para clustering
    X        = X,         # matriz de atributos
    D0e      = D0e,       # distancia socioeconómica y estructural
    D1e      = D1e        # distancia espacial
  )
}

# Ejecución de la preparación y fusión de matrices
message("\n--- PREPARACIÓN Y FUSIÓN DE MATRICES ---")
res_final      <- preparar_y_fusionar(df0, w_espacial = W_ESP)
Dfused_e_final <- res_final$Dfused_e
X_final        <- res_final$X
D0e_final      <- res_final$D0e
D1e_final      <- res_final$D1e
# -------------------------------------------------------------------
# Visualización de la matriz fusionada y exportación de las matrices
# -------------------------------------------------------------------
p_dist_fused <- factoextra::fviz_dist(
  Dfused_e_final, order = TRUE, show_labels = TRUE, lab_size = 2,
  gradient = list(low="blue", mid="white", high="red")
) + ggplot2::ggtitle(sprintf("Dfused_e – Fusión Final (w_atr=%.2f, w_esp=%.2f)", W_ATR, W_ESP))
print(p_dist_fused)

write_dist_csv(Dfused_e_final, file.path(OUT_DIR, "Dfused_e_fusion_final.csv"))
write_dist_csv(D0e_final,      file.path(OUT_DIR, "D0e_gower_atributos.csv"))
write_dist_csv(D1e_final,      file.path(OUT_DIR, "D1e_espacial_euclidiana.csv"))
# ---------------------------------------
# Exploración del número de k-clústeres 
# ----------------------------------------
stopifnot(inherits(Dfused_e_final, "dist"))
nobj <- attr(Dfused_e_final, "Size")
KS   <- 2:min(K_MAX, nobj - 1)
K_FINAL_MAX <- max(KS)

# Dendrograma Ward
hcF <- hclust(Dfused_e_final, method = "ward.D2")
print_quiet(
  factoextra::fviz_dend(hcF, k = NULL, rect = FALSE, cex = 0.6) +
    ggplot2::ggtitle("Dendrograma Ward.D2 – Exploración sobre Dfused_e_final")
)

# Gráficos (codo / silhouette)
# Codo
print_quiet(
  factoextra::fviz_nbclust(as.matrix(Dfused_e_final), FUNcluster = cluster::pam,
                           method = "wss", k.max = K_FINAL_MAX, diss = Dfused_e_final) +
    ggplot2::ggtitle("Elbow (WSS) – PAM sobre Matriz Fusionada")
)
# Silhouette
print_quiet(
  factoextra::fviz_nbclust(as.matrix(Dfused_e_final), FUNcluster = cluster::pam,
                           method = "silhouette", k.max = K_FINAL_MAX, diss = Dfused_e_final) +
    ggplot2::ggtitle("Silhouette – PAM sobre Matriz Fusionada")
)
# valores de codo
elbow_vals <- data.frame(k = KS, wss_objetivo = NA_real_)
for (i in seq_along(KS)) {
  k_try   <- KS[i]
  pam_res <- cluster::pam(Dfused_e_final, k = k_try, diss = TRUE)
  elbow_vals$wss_objetivo[i] <- pam_res$objective[1]
}
cat("\n--- VALORES DEL CODO (WSS/objetivo PAM) POR k ---\n")
print(elbow_vals)
# Indice de silhouette
sil_vals <- data.frame(k = KS, silhouette_avg = NA_real_)
for (i in seq_along(KS)) {
  k_try   <- KS[i]
  pam_res <- cluster::pam(Dfused_e_final, k = k_try, diss = TRUE)
  sil_res <- cluster::silhouette(pam_res$clustering, Dfused_e_final)
  sil_vals$silhouette_avg[i] <- mean(sil_res[, 3])
}
cat("\n--- VALORES DE SILHOUETTE PROMEDIO POR k ---\n")
print(sil_vals)
#
top3_k <- sil_vals %>%
  dplyr::arrange(dplyr::desc(silhouette_avg)) %>%
  dplyr::slice_head(n = 3)

cat("\n--- TOP 3 VALORES DE k SEGÚN SILHOUETTE ---\n")
print(top3_k)

readr::write_csv(top3_k,  file.path(OUT_DIR, "top3_k_silueta.csv"))

message("La selección final de k para análisis detallado se controla con K_ELEGIDOS.")
# --------------------------
# Bucle PAM para K_ELEGIDOS
# --------------------------
K_ELEGIDOS <- sort(unique(K_ELEGIDOS))
K_ELEGIDOS <- K_ELEGIDOS[K_ELEGIDOS >= 2 & K_ELEGIDOS <= min(K_MAX, nobj - 1)]

if (length(K_ELEGIDOS) == 0) stop("K_ELEGIDOS no contiene valores válidos.")
message(sprintf("\nSe analizarán en detalle los k: %s", paste(K_ELEGIDOS, collapse = ", ")))

for (k_cur in K_ELEGIDOS) {

  message("\n===================================================")
  message(sprintf("========== PAM + RESÚMENES PARA k = %d ==========", k_cur))
  message("===================================================\n")

  # (1) Se ejecuta PAM
  pam_k <- cluster::pam(Dfused_e_final, k = k_cur, diss = TRUE)

  # (2) Exportación: asignación de barrios
  asignacion_k <- data.frame(
    BARRIO  = rownames(df0),
    cluster = pam_k$clustering
  )
  fname_asig <- sprintf("asignacion_barrios_k%d.csv", k_cur)
  fpath_asig <- file.path(OUT_DIR, fname_asig)
  utils::write.csv(asignacion_k, fpath_asig, row.names = FALSE)

  message(sprintf("[OK] Export asignación (original): %s", file.path(getwd(), fpath_asig)))

  # (3) Resumen por clúster
  df_resumen_k <- df0
  df_resumen_k$cluster_id <- factor(pam_k$clustering)

  resumen_k <- df_resumen_k %>%
    dplyr::group_by(cluster_id) %>%
    dplyr::summarise(
      n_barrios        = dplyr::n(),
      poblacion_total  = sum(.data[[map$poblacion]], na.rm = TRUE),
      ipm_media        = mean(.data[[map$ipm]],       na.rm = TRUE),
      danio_media      = mean(.data[[map$danio]],     na.rm = TRUE),
      estrato_moda     = moda_numerica(.data[[map$estrato]]),
      demanda_total    = sum(DEMANDA, na.rm = TRUE),    
      .groups = "drop"
    ) %>%
    dplyr::arrange(cluster_id)

  message(sprintf("\n--- Resumen por clúster (k = %d) ---", k_cur))
  print(resumen_k)
  # Exportar resumen por clúster
  readr::write_csv(resumen_k, file.path(OUT_DIR, sprintf("resumen_cluster_k%d.csv", k_cur)))

  # (4) Perfil normalizado (GRÁFICO) 
  DANIO_MIN_G <- min(as.numeric(df0[[map$danio]]),     na.rm = TRUE)
  DANIO_MAX_G <- max(as.numeric(df0[[map$danio]]),     na.rm = TRUE)
  
  IPM_MIN_G   <- min(as.numeric(df0[[map$ipm]]),       na.rm = TRUE)
  IPM_MAX_G   <- max(as.numeric(df0[[map$ipm]]),       na.rm = TRUE)
  
  POB_MIN_G   <- min(as.numeric(df0[[map$poblacion]]), na.rm = TRUE)
  POB_MAX_G   <- max(as.numeric(df0[[map$poblacion]]), na.rm = TRUE)
  
  EST_MIN_G <- 1
  EST_MAX_G <- 6
  
  perfil_k <- df_resumen_k %>%
    dplyr::group_by(cluster_id) %>%
    dplyr::summarise(
      danio_mean     = mean(.data[[map$danio]],     na.rm = TRUE),
      ipm_mean       = mean(.data[[map$ipm]],       na.rm = TRUE),
      poblacion_mean = mean(.data[[map$poblacion]], na.rm = TRUE),
      estrato_moda   = moda_numerica(.data[[map$estrato]]),
      .groups = "drop"
    ) %>%
    dplyr::mutate(
      estrato_num = as.numeric(sub("/.*", "", estrato_moda)),
      estrato_num = pmin(pmax(estrato_num, EST_MIN_G), EST_MAX_G)      
      )

  perfil_long <- perfil_k %>%
    dplyr::mutate(
      danio     = norm01_lim(danio_mean,     DANIO_MIN_G, DANIO_MAX_G),
      ipm       = norm01_lim(ipm_mean,       IPM_MIN_G,   IPM_MAX_G),
      poblacion = norm01_lim(poblacion_mean, POB_MIN_G,   POB_MAX_G),
      estrato   = 1 - norm01_lim(estrato_num, EST_MIN_G, EST_MAX_G)
    ) %>%
    dplyr::select(cluster_id, danio, ipm, estrato, poblacion) %>%
    tidyr::pivot_longer(cols = -cluster_id, names_to = "variable", values_to = "valor")

  perfil_long$variable <- factor(
    perfil_long$variable,
    levels = c("danio","ipm","estrato","poblacion"),
    labels = c("% DAÑO","IPM","ESTRATO (invertido)","POBLACIÓN")
  )

  print(
    ggplot2::ggplot(perfil_long,
                    ggplot2::aes(x = variable, y = valor, group = cluster_id, colour = cluster_id)) +
      ggplot2::geom_line() +
      ggplot2::geom_point(size = 2) +
      ggplot2::scale_y_continuous(limits = c(0,1)) +
      ggplot2::labs(
        title  = paste0("Perfil normalizado de clústeres – k = ", k_cur),
        x = "Variable",
        y = "Valor normalizado (0–1)",
        colour = "Clúster"
      ) +
      ggplot2::theme_minimal()
  )
}
# ------------------------------------------------------------
# Datos base para cálculo de radio y cobertura (por clúster)
# ------------------------------------------------------------
# Se construye una tabla base por barrio con variables que afectan la
# vulnerabilidad vial y, por ende, el radio efectivo de atención:
#  - ESTRATO: condiciones del entorno (vulnerabilidad relativa)
#  - DANO: % de daño estimado (impacto esperado en infraestructura/vías)
# Estas variables se agregan a nivel de clúster para calcular un Vi_clúster.
cov_base <- data.frame(
  BARRIO  = rownames(df0),
  ESTRATO = as.integer(df0[[map$estrato]]),
  DANO    = as.numeric(df0[[map$danio]]),
  stringsAsFactors = FALSE
)

# Valores globales (mín/máx) del % de daño para normalización min–max a nivel ciudad
# (se usan luego para construir V_DANO en el cálculo del radio por clúster).
DANO_MIN <- min(cov_base$DANO, na.rm = TRUE)
DANO_MAX <- max(cov_base$DANO, na.rm = TRUE)
# ------------------------------------------------------------
# Bucle por escenario k seleccionado (K_ELEGIDOS)
# ------------------------------------------------------------
# Para cada k:
#  1) Lee la asignación barrio→clúster (output de la fase PAM)
#  2) Localiza un CMT por clúster usando centroide ponderado por DEMANDA
#  3) Calcula radio de cobertura por clúster con T_RTA y V_NOM ajustados por vulnerabilidad
#  4) Exporta salidas requeridas
for (k_cur in K_ELEGIDOS) {
  
  # Archivo esperado con la asignación de clústeres para este k
  asig_path <- file.path(OUT_DIR, sprintf("asignacion_barrios_k%d.csv", k_cur))
  if (!file.exists(asig_path)) {
    warning(sprintf("[WARN] No existe %s. Salto k=%d.", asig_path, k_cur))
    next
  }
  
  # Lectura de asignación barrio-clúster
  asig <- utils::read.csv(asig_path, stringsAsFactors = FALSE)
  
  # Validación mínima de estructura: se requieren columnas BARRIO y cluster
  stopifnot(all(c("BARRIO","cluster") %in% names(asig)))
  
  # ----------------------------------------------------------
  # (1) Localización de CMT por clúster: centroide ponderado por DEMANDA
  # ----------------------------------------------------------
  # Se integra la asignación del clúster con la tabla barrios_dem (que contiene
  # DEMANDA y coordenadas X,Y por barrio)
  df_loc <- asig %>%
    dplyr::left_join(barrios_dem, by = "BARRIO") %>%
    dplyr::filter(!is.na(DEMANDA), !is.na(X), !is.na(Y))
  
  # Para cada clúster, el CMT se ubica en el centroide ponderado:
  centros <- df_loc %>%
    dplyr::group_by(cluster) %>%
    dplyr::summarise(
      total_demanda = sum(DEMANDA, na.rm = TRUE),
      lon_centroide = sum(X * DEMANDA, na.rm = TRUE) / total_demanda,
      lat_centroide = sum(Y * DEMANDA, na.rm = TRUE) / total_demanda,
      .groups = "drop"
    )
  
  # Exportación requerida (salida original): centros por clúster
  out_centros <- file.path(OUT_DIR, sprintf("centros_medicos_centroides_k%d.csv", k_cur))
  utils::write.csv(centros, out_centros, row.names = FALSE)
  message(sprintf("[OK] k=%d -> exportado (original): %s", k_cur, file.path(getwd(), out_centros)))
  # ----------------------------------------------------------
  # (2) Cálculo del radio por clúster (para evaluación de cobertura)
  # ----------------------------------------------------------
  # Se arma una base barrio-clúster y se adjuntan ESTRATO y DANO por barrio.
  # Esto permite agregar a nivel de clúster: estrato representativo (moda)
  # y daño promedio, que se usan como aproximación de vulnerabilidad del clúster.
  df_cov <- asig %>%
    dplyr::transmute(
      BARRIO  = as.character(BARRIO),
      CLUSTER = as.integer(cluster)
    ) %>%
    dplyr::left_join(cov_base, by = "BARRIO")
  
  # Resumen por clúster con parámetros base del modelo (T_RTA y V_NOM) y agregados:
  resumen_cob <- df_cov %>%
    dplyr::group_by(CLUSTER) %>%
    dplyr::summarise(
      T_RTA_MIN = T_RTA_MIN,
      V_NOM_KMH = V_NOM_KMH,
      ESTRATO   = moda_estrato(ESTRATO),
      DANO_PCT  = mean(DANO, na.rm = TRUE),
      .groups   = "drop"
    ) %>%
    dplyr::mutate(
      # Normalización de vulnerabilidad (0–1) para estrato y daño
      V_ESTRATO = (6 - ESTRATO) / 5,
      V_DANO    = (DANO_PCT - DANO_MIN) / (DANO_MAX - DANO_MIN),
      
      # Índice de vulnerabilidad del clúster (ejemplo: ponderación 50/50)
      Vi        = 0.5 * (V_ESTRATO + V_DANO),
      
      # Velocidad efectiva ajustada por vulnerabilidad (reduce V_NOM)
      Vf_kmh    = V_NOM_KMH * (1 - Vi),
      
      # Radio efectivo (m): distancia = velocidad (km/h) * tiempo (h) -> km -> m
      RADIO_M   = Vf_kmh * (T_RTA_MIN / 60) * 1000
    ) %>%
    # Se adjuntan las coordenadas del centroide del CMT por clúster
    dplyr::left_join(
      centros %>%
        dplyr::transmute(
          CLUSTER = as.integer(cluster),
          lon_centroide, lat_centroide, total_demanda
        ),
      by = "CLUSTER"
    ) %>%
    dplyr::arrange(CLUSTER)
  # ---------------------------
  # (3) Exportación de radios 
  # --------------------------

  out_radios <- file.path(OUT_DIR, sprintf("radio_cluster_k%d.csv", k_cur))
  utils::write.csv(resumen_cob, out_radios, row.names = FALSE)
  
  message(sprintf("[OK] k=%d -> exportado (original): %s", k_cur, file.path(getwd(), out_radios)))
}
# ------------------------------------
# COBERTURA FINAL (DEMANDA CUBIERTA) 
# ------------------------------------
# CRS de entrada (coordenadas X,Y del consolidado) y CRS métrico para distancias
CRS_INPUT  <- CRS_GEO
CRS_METRIC <- CRS_PROJ

# ------------------------------------------------------------
# Función: calcula cobertura para un escenario k
# ------------------------------------------------------------
# Para un k dado:
#  1) Lee asignación barrio→clúster, centros por clúster y radio por clúster
#  2) Calcula distancia barrio→centro (en metros)
#  3) Determina cobertura: CUBIERTO = (DIST_M <= RADIO_M)
#  4) Exporta resumen (por clúster + total ciudad)
cobertura_por_k <- function(k_cur){
  
  # Rutas de archivos requeridos
  asig_path    <- file.path(OUT_DIR, sprintf("asignacion_barrios_k%d.csv", k_cur))
  centros_path <- file.path(OUT_DIR, sprintf("centros_medicos_centroides_k%d.csv", k_cur))
  radios_path  <- file.path(OUT_DIR, sprintf("radio_cluster_k%d.csv", k_cur))
  
  # Lectura de inputs del escenario
  asig    <- utils::read.csv(asig_path, stringsAsFactors = FALSE)
  centros <- utils::read.csv(centros_path, stringsAsFactors = FALSE)
  radios  <- utils::read.csv(radios_path, stringsAsFactors = FALSE)
 
  # Estandarización de campos:
    # Se homogeneizan nombres y tipos para asegurar uniones consistentes.
  asig <- asig %>%
    dplyr::transmute(
      BARRIO  = as.character(BARRIO),
      CLUSTER = as.integer(cluster)
    )

  centros <- centros %>%
    dplyr::transmute(
      CLUSTER = as.integer(cluster),
      lon_centroide = as.numeric(lon_centroide),
      lat_centroide = as.numeric(lat_centroide),
      total_demanda_centro = as.numeric(total_demanda)
    )

  radios <- radios %>%
    dplyr::transmute(
      CLUSTER = as.integer(CLUSTER),
      RADIO_M = as.numeric(RADIO_M)
    )
  # --------------------------------------
  # Construcción de base barrio-CMT-radio
  # --------------------------------------
  # Se une:
  #  - Asignación (barrio-clúster)
  #  - Barrios con DEMANDA y coordenadas X,Y (barrios_dem)
  #  - Centroide por clúster
  #  - Radio por clúster
  
  base <- asig %>%
    dplyr::left_join(barrios_dem, by = "BARRIO") %>%
    dplyr::filter(!is.na(CLUSTER), !is.na(DEMANDA), !is.na(X), !is.na(Y)) %>%
    dplyr::left_join(centros, by = "CLUSTER") %>%
    dplyr::left_join(radios,  by = "CLUSTER") %>%
    dplyr::filter(!is.na(lon_centroide), !is.na(lat_centroide), !is.na(RADIO_M))
  
  
  #Distancia barrio-centro (en metros)
    # Se convierten barrios y centros a objetos sf y se reproyectan a CRS métrico
    # para que la distancia sea calculada en unidades métricas (m).
  sf_barrios <- sf::st_as_sf(base, coords = c("X","Y"), crs = CRS_INPUT, remove = FALSE) %>%
    sf::st_transform(CRS_METRIC)

  sf_centros <- sf::st_as_sf(base, coords = c("lon_centroide","lat_centroide"), crs = CRS_INPUT, remove = FALSE) %>%
    sf::st_transform(CRS_METRIC)
  
  # Distancia por elemento (cada barrio contra su propio centro de clúster)
  dist_m <- as.numeric(sf::st_distance(sf_barrios, sf_centros, by_element = TRUE))

  # --------------------------
  # Cobertura a nivel barrio
  # --------------------------
  detalle <- base %>%
    dplyr::mutate(
      DIST_M   = dist_m,
      CUBIERTO = DIST_M <= RADIO_M
    )
  # Resumen por clúster
  resumen_cluster <- detalle %>%
    dplyr::group_by(CLUSTER) %>%
    dplyr::summarise(
      ESCENARIO           = paste0("k", k_cur),
      N_BARRIOS_CLUSTER   = dplyr::n(),
      N_BARRIOS_CUBIERTOS = sum(CUBIERTO),
      DEMANDA_TOTAL       = sum(DEMANDA, na.rm = TRUE),
      DEMANDA_CUBIERTA    = sum(DEMANDA[CUBIERTO], na.rm = TRUE),
      PROP_COBERTURA      = ifelse(DEMANDA_TOTAL > 0, DEMANDA_CUBIERTA / DEMANDA_TOTAL, NA_real_),
      .groups = "drop"
    ) %>%
    dplyr::arrange(CLUSTER)
  
  # Resumen total ciudad 
  resumen_ciudad <- resumen_cluster %>%
    dplyr::summarise(
      ESCENARIO           = paste0("k", k_cur),
      CLUSTER             = "TOTAL_CIUDAD",
      N_BARRIOS_CLUSTER   = sum(N_BARRIOS_CLUSTER),
      N_BARRIOS_CUBIERTOS = sum(N_BARRIOS_CUBIERTOS),
      DEMANDA_TOTAL       = sum(DEMANDA_TOTAL),
      DEMANDA_CUBIERTA    = sum(DEMANDA_CUBIERTA),
      PROP_COBERTURA      = ifelse(DEMANDA_TOTAL > 0, DEMANDA_CUBIERTA / DEMANDA_TOTAL, NA_real_)
    )

  # Consolidación final (clústeres + total ciudad)
  resumen_final <- dplyr::bind_rows(
    resumen_cluster %>% dplyr::mutate(CLUSTER = as.character(CLUSTER)),
    resumen_ciudad
  )

  #Exportar resúmen creado
  readr::write_csv(resumen_final, file.path(OUT_DIR, sprintf("cobertura_resumen_k%d.csv", k_cur)))
  message(sprintf("[OK] Cobertura calculada y exportada para k=%d", k_cur))
}
# Ejecución de cobertura para cada escenario k seleccionado
for (k_cur in K_ELEGIDOS) {
  cobertura_por_k(k_cur)
}

message("\n[FIN] Modelo ejecutado: clusterización + localización + radio + cobertura (con exports completos).")
